SECTION 1: DATA CURATION

1.1 Data Files

#setwd('Dropbox/Mycobiome Xav Analysis/')

######### LIBRARIES 

        library(phyloseq)
        library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.1.0     ✓ dplyr   1.0.5
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
        library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-7
        library(brms)
## Loading required package: Rcpp
## Loading 'brms' package (version 2.15.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
## 
## Attaching package: 'brms'
## The following object is masked from 'package:phyloseq':
## 
##     nsamples
## The following object is masked from 'package:stats':
## 
##     ar
        library(dplyr)
        library(tidybayes)
## 
## Attaching package: 'tidybayes'
## The following objects are masked from 'package:brms':
## 
##     dstudent_t, pstudent_t, qstudent_t, rstudent_t
        library(tidyr)
        library(cowplot)
        library(microbiome)
## 
## microbiome R package (microbiome.github.com)
##     
## 
## 
##  Copyright (C) 2011-2020 Leo Lahti, 
##     Sudarshan Shetty et al. <microbiome.github.io>
## 
## Attaching package: 'microbiome'
## The following object is masked from 'package:vegan':
## 
##     diversity
## The following object is masked from 'package:ggplot2':
## 
##     alpha
## The following object is masked from 'package:base':
## 
##     transform
        options(mc.cores=parallel::detectCores())
        

#Load Data Files
    myco<-readRDS('clean.its.rds')
    micro<-readRDS('clean.16s.rds')
    
##Metadata
    metadata<-read.csv('Sample Metadata with Diet May20.csv',header=T)
        #head(metadata)   
    
##RNG SEED
    set.seed(21042020)
    
## Use multiple Cores for Model Fitting (brms)    
       options(mc.cores=parallel::detectCores()) 
    
       
#Global Plot Options
            plotopts<- theme(axis.text=element_text(size=12),axis.title=element_text(size=16),strip.text=element_text(size=20),legend.title = element_text(size=18),legend.text = element_text(size=18))

SECTION 2 NETWORK ANALYSIS

2.1 Correlations Among Microbial ASVs

#library(devtools)
#install_github("zdk123/SpiecEasi")


### Core Datasets of Myco and Micro Sample 
microbial_common_samples<-Reduce(intersect, list(rownames(sample_data(myco)),rownames(sample_data(micro))))

myco_common<-subset_samples(myco,Sample %in% microbial_common_samples)
micro_common<-subset_samples(micro,Sample %in% microbial_common_samples)

#Checks 
phyloseq::nsamples(myco_common) == phyloseq::nsamples(micro_common)
## [1] TRUE
all.equal(as.character(sample_data(myco_common)$Sample),as.character(sample_data(micro_common)$Sample)) #TRUE
## [1] TRUE
#######  MERGING OBJECTS      
#rename if required and check again
taxa_names(micro_common) <- paste0("bacteria-", taxa_names(micro_common))
#taxa_names(micro_common)

taxa_names(myco_common) <- paste0("fungi-", taxa_names(myco_common))
#taxa_names(myco_common)

#Make Sure No Duplicate SV Names
  table(taxa_names(micro_common) %in% taxa_names(myco_common)) ## ALL FALSE
## 
## FALSE 
##   803
#Merge Objects
  microbial_merge<-merge_phyloseq(micro_common,myco_common)

#More Checking - Samples and Microbial Taxa Totals
ntaxa(microbial_merge) == (ntaxa(micro_common) + ntaxa(myco_common)) #TRUE
## [1] TRUE
phyloseq::nsamples(microbial_merge) == length(microbial_common_samples) #TRUE
## [1] TRUE
### Metadata Checking
table(sample_data(microbial_merge)$Common_name)
## 
##  African palm weevil larvae                   Bank vole 
##                           6                          10 
##               Barred hamlet                    Blackcap 
##                           5                           5 
##                    Blue tit                Brown shrimp 
##                           3                           2 
##             Capuchin monkey                Carrion crow 
##                           9                           2 
##                  Chiffchaff                   Cockroach 
##                           3                           7 
##               Collared dove         Common midwife toad 
##                           5                          11 
##                 Common toad         Eastern black rhino 
##                           2                          10 
##                European eel       Foureye butterflyfish 
##                           1                           9 
##                   Goldfinch                   Great tit 
##                           3                           3 
## Greater white-toothed shrew               Grey squirrel 
##                           7                           5 
##                   Hard tick                    Hedgehog 
##                           2                          11 
##                   Honey bee        Lesser horseshoe bat 
##                           6                           2 
##   Light-bellied brent goose            Northern muriqui 
##                           9                          10 
##          Phofung river frog                 Pygmy shrew 
##                           4                          10 
##                    Red deer                Red squirrel 
##                          10                          11 
##                Reed bunting                Reed warbler 
##                           2                           2 
##                       Robin                    Roe deer 
##                           4                           4 
##                 Song thrush                  Stock dove 
##                           5                           1 
##         Striped field mouse                  Tsetse fly 
##                          10                           6 
##                 Turtle dove                 Vase sponge 
##                           7                          10 
##                   Wild pony                  Wood mouse 
##                          10                           9 
##                  Woodpigeon         Yellow-necked mouse 
##                           4                          10 
##                Yellowhammer           Yellowhead wrasse 
##                           5                           5

2.2 Fitting Co-Occurrence Models

########## SPIEC-EASI MODELSs         

library(SpiecEasi)
library(igraph)

### Rename Phyla
tax_table(microbial_merge)[,2][which(tax_table(microbial_merge)[,2]=="Kingdom_k__Fungi")] <-"Kingdom_Fungi"
tax_table(microbial_merge)[,2]<-gsub("^p__","",tax_table(microbial_merge)[,2])


## Function for Colour Edges by Sign of Relationship 
extract_edge_cols<-function(physeq_obj,speic_obj,igraph_obj){
  
  otu.ids=rownames(tax_table(physeq_obj))
  betaMat=as.matrix(symBeta(getOptBeta(speic_obj)))
  
  edges=E(igraph_obj)
  edge.colors=c()
  for(e.index in 1:length(edges)){
    adj.nodes=ends(igraph_obj,edges[e.index])
    xindex=which(otu.ids==adj.nodes[1])
    yindex=which(otu.ids==adj.nodes[2])
    beta=betaMat[xindex,yindex]
    if(beta>0){
      edge.colors=append(edge.colors,"forestgreen")
    }else if(beta<0){
      edge.colors=append(edge.colors,"red")
    }
  }
  
  return (edge.colors)
}



### Palette for Kingdom 
library(RColorBrewer)   
length(unique(tax_table(microbial_merge)[,2]))  #21 Phyla
## [1] 21
global_phy_cols<-data.frame(Phylum=unique(tax_table(microbial_merge)[,2]))
#global_phy_cols$color<-colorRampPalette(brewer.pal(11,"Spectral"))(nrow(global_phy_cols))
global_phy_cols$color<-c(brewer.pal(12,"Paired"),brewer.pal(9,"Pastel1"))



### Mammal Test 
  mammals <- subset_samples(microbial_merge,Taxa=="Mammal") %>% filter_taxa(function(x) sum(x > 3) > (0.05*length(x)), TRUE) 
  mammal_speic<- mammals %>% spiec.easi(method='mb', lambda.min.ratio=1e-2, nlambda=20,pulsar.params=list(rep.num=50))
## Applying data transformations...
## Selecting model with pulsar using stars...
## Fitting final estimate with mb...
## done
  ig.mammal <- adj2igraph(getRefit(mammal_speic),vertex.attr=list(name=taxa_names(mammals),Kingdom=tax_table(mammals)[,1],Phylum=tax_table(mammals)[,2]))

#Assign Shapes to Kingdoms + Colours to Phylum
V(ig.mammal)$shape<-ifelse(V(ig.mammal)$Kingdom=="Bacteria","circle","square")
V(ig.mammal)$color<-global_phy_cols$color[match(V(ig.mammal)$Phylum,global_phy_cols$Phylum)]

### Assign Edge Colours Based on Positive Or Negativ         
mammal_edgecols<-extract_edge_cols(mammals,mammal_speic,ig.mammal)          
E(ig.mammal)$color=mammal_edgecols   
#E(ig.mammal)$weight<-3

#Plot   
plot(ig.mammal,vertex.size=9,vertex.label=NA)

plot_network(ig.mammal)

### Birds 
birds <- subset_samples(microbial_merge,Taxa=="Bird") %>% filter_taxa(function(x) sum(x > 3) > (0.05*length(x)), TRUE) 
bird_speic<- birds %>% spiec.easi(method='mb', lambda.min.ratio=1e-2, nlambda=20,pulsar.params=list(rep.num=50))
## Applying data transformations...
## Selecting model with pulsar using stars...
## Fitting final estimate with mb...
## done
ig.bird <- adj2igraph(getRefit(bird_speic),vertex.attr=list(name=taxa_names(birds),Kingdom=tax_table(birds)[,1],Phylum=tax_table(birds)[,2]))


#Assign Shapes to Kingdoms + Colours to Phylum
V(ig.bird)$shape<-ifelse(V(ig.bird)$Kingdom=="Bacteria","circle","square")
V(ig.bird)$color<-global_phy_cols$color[match(V(ig.bird)$Phylum,global_phy_cols$Phylum)]

### Assign Edge Colours Based on Positive Or Negativ         
bird_edgecols<-extract_edge_cols(birds,bird_speic,ig.bird)          
E(ig.bird)$color=bird_edgecols   

#Plot
plot(ig.bird,vertex.size=9,vertex.label=NA,circular=T)

### Amphibians 
amphibians <- subset_samples(microbial_merge,Taxa=="Amphibian") %>% filter_taxa(function(x) sum(x > 3) > (0.05*length(x)), TRUE) 
amphib_speic<- amphibians %>% spiec.easi(method='mb', lambda.min.ratio=1e-2, nlambda=20,pulsar.params=list(rep.num=50))
## Applying data transformations...
## Selecting model with pulsar using stars...
## Fitting final estimate with mb...
## done
ig.amphib <- adj2igraph(getRefit(amphib_speic),vertex.attr=list(name=taxa_names(amphibians),Kingdom=tax_table(amphibians)[,1],Phylum=tax_table(amphibians)[,2]))


#Assign Shapes to Kingdoms + Colours to Phylum
V(ig.amphib)$shape<-ifelse(V(ig.amphib)$Kingdom=="Bacteria","circle","square")
V(ig.amphib)$color<-global_phy_cols$color[match(V(ig.amphib)$Phylum,global_phy_cols$Phylum)]

### Assign Edge Colours Based on Positive Or Negativ         
amphibian_edgecols<-extract_edge_cols(amphibians,amphib_speic,ig.amphib)          
E(ig.amphib)$color=amphibian_edgecols   

#Plot
plot(ig.amphib,vertex.size=9,vertex.label=NA,main="Amphibians",circular=T)

#legend("topleft",legend=phylumcols_amphib[,1],fill=phylumcols_amphib[,2])

#### Fish 
fish <- subset_samples(microbial_merge,Taxa=="Fish") %>% filter_taxa(function(x) sum(x > 3) > (0.05*length(x)), TRUE) 
fish_speic<- fish %>% spiec.easi(method='mb', lambda.min.ratio=1e-2, nlambda=20,pulsar.params=list(rep.num=50))
## Applying data transformations...
## Selecting model with pulsar using stars...
## Fitting final estimate with mb...
## done
ig.fish <- adj2igraph(getRefit(fish_speic),vertex.attr=list(name=taxa_names(fish),Kingdom=tax_table(fish)[,1],Phylum=tax_table(fish)[,2]))


#Assign Shapes to Kingdoms + Colours to Phylum
V(ig.fish)$shape<-ifelse(V(ig.fish)$Kingdom=="Bacteria","circle","square")
V(ig.fish)$color<-global_phy_cols$color[match(V(ig.fish)$Phylum,global_phy_cols$Phylum)]


### Assign Edge Colours Based on Positive Or Negativ         
fish_edgecols<-extract_edge_cols(fish,fish_speic,ig.fish)          
E(ig.fish)$color=fish_edgecols            

#Plot
plot(ig.fish,vertex.size=9,vertex.label=NA,main="Fish",circular=T)

#legend("topleft",legend=phylumcols_amphib[,1],fill=phylumcols_amphib[,2])

#### Insecta 
sample_data(microbial_merge)$Class<-metadata$Class[match(sample_data(microbial_merge)$Common_name,metadata$Common_name)]      
inverts <- subset_samples(microbial_merge,Class=="Insecta") %>% filter_taxa(function(x) sum(x > 3) > (0.05*length(x)), TRUE) 
invert_speic<- inverts %>% spiec.easi(method='mb', lambda.min.ratio=1e-2, nlambda=20,pulsar.params=list(rep.num=50))
## Applying data transformations...
## Selecting model with pulsar using stars...
## Fitting final estimate with mb...
## done
ig.invert <- adj2igraph(getRefit(invert_speic),vertex.attr=list(name=taxa_names(inverts),Kingdom=tax_table(inverts)[,1],Phylum=tax_table(inverts)[,2]))


#Assign Shapes to Kingdoms + Colours to Phylum
V(ig.invert)$shape<-ifelse(V(ig.invert)$Kingdom=="Bacteria","circle","square")
V(ig.invert)$color<-global_phy_cols$color[match(V(ig.invert)$Phylum,global_phy_cols$Phylum)]

### Assign Edge Colours Based on Positive Or Negativ         
invert_edgecols<-extract_edge_cols(inverts,invert_speic,ig.invert)          
E(ig.invert)$color=invert_edgecols   
#Plot
plot(ig.invert,vertex.size=9,vertex.label=NA,main="Insecta")

plot(ig.invert,vertex.size=9,vertex.label=NA,main="Invertebrate",circular=TRUE)

2.3 Plot of All Networks

         # jpeg(filename="network_plots.jpeg", width=300, height=320, units="mm", res=600)
            final.fig <- par(mfcol=c(3, 2),mar=c(1,1,1,1),oma=c(1,1,1,1))
            plot(ig.mammal,vertex.size=12,vertex.label=NA,main="Mammalia",circular=T,label.cex=10)
            plot(ig.bird,vertex.size=12,vertex.label=NA,circular=T,main="Aves",circular=T)
            plot(ig.amphib,vertex.size=12,vertex.label=NA,main="Amphibia",circular=T)
            plot(ig.fish,vertex.size=12,vertex.label=NA,main="Actinopterygii",circular=T)
            plot(ig.invert,vertex.size=12,vertex.label=NA,main="Insecta",circular=TRUE)
            plot.new()
            legend("topleft",legend=global_phy_cols[,1],fill=global_phy_cols[,2],border=NULL,cex=1.8,horiz=FALSE,ncol=2)
            legend("bottomleft",legend=c("Bacteria","Fungi"),pch=c(21,22),cex=2.5,pt.cex=3)

          #   par(final.fig)
          # dev.off()

2.4 Summary Stats of Network Interactions

#Assemble Data
      taxa_edge_data<-data.frame(cols=c(mammal_edgecols,bird_edgecols,amphibian_edgecols,fish_edgecols,invert_edgecols),taxa=rep(c("Mammal","Bird","Amphibian","Fish","Invert"),times=c(length(mammal_edgecols),length(bird_edgecols),length(amphibian_edgecols),length(fish_edgecols),length(invert_edgecols))))

taxa_edge_data<- taxa_edge_data %>% mutate(sign=ifelse(cols=="forestgreen","Positive","Negative"))

taxa_edge_data %>% group_by(taxa) %>% summarise(npos=length(which(sign=="Positive")),nneg=length(which(sign=="Negative")))
## # A tibble: 5 x 3
##   taxa       npos  nneg
##   <chr>     <int> <int>
## 1 Amphibian    68    10
## 2 Bird        241     6
## 3 Fish         64     1
## 4 Invert      121     0
## 5 Mammal      686     9

2.5 Co-Occurrence Stats

2.5.1 Functions for Bacteria-Fungus Co-Occurrence Analysis

library(igraph)

#FUnction to Pull Out Edges and Filter to Bacteria-Fungi Interactions 
    edge_extract<-function(igraphobject){
      
    #Pull Out Edge and Vertex Data  
      edges_data_frame <- get.data.frame(igraphobject, what = "edges")
      v_data_frame <- get.data.frame(igraphobject, what = "vertices")
    
    #Annotate With Metadata    
      edges_data_frame$king1<-v_data_frame$Kingdom[match(edges_data_frame$from,v_data_frame$name)]
      edges_data_frame$king2<-v_data_frame$Kingdom[match(edges_data_frame$to,v_data_frame$name)]
        
    edges_data_frame$kingpaste<-with(edges_data_frame,paste(king1,king2))
    edges_data_frame$phylum1<-v_data_frame$Phylum[match(edges_data_frame$from,v_data_frame$name)]
        edges_data_frame$phylum2<-v_data_frame$Phylum[match(edges_data_frame$to,v_data_frame$name)]
  #Filter to  Postive Co-Occurence        
    edges_bf<-edges_data_frame %>% filter(kingpaste=="Bacteria k__Fungi" & color=="forestgreen") %>% filter(phylum2!="Kingdom_Fungi")
    return(edges_bf)
    }


#### FUnction to Tabulate and Provide a Base Heatmap 
  heatplot<-function(edge_extract_output){
    
  #Get Data   
    d1<-with(edge_extract_output,table(phylum1,phylum2)) %>% as.data.frame() %>% mutate(Prop=Freq/sum(Freq))
    p1<- d1 %>% ggplot(aes(phylum1,phylum2,fill=Prop)) + geom_tile() + labs(x="Bacterial Phylum",y="Fungal Phylum") + plotopts + theme(axis.text.x = element_text(angle=45,hjust = 1))
    return(p1)
  }  
  
    

#### Function to Permute Frequency of Co-occurrences   

  occur_rand<-function(edge_extract_output,nsims){
    
    #Work Out Which is The Most Common 
      d1<-with(edge_extract_output,table(phylum1,phylum2)) %>% as.data.frame() %>% arrange(desc(Freq))
      
    #Permutations
      #Empty Vecgtor for Frequency 
        permute_vec<-numeric(nsims)
          for(k in 1:nsims){
            #Shuffle Phylum 2
              edge_extract_output$phylum2shuffle<-sample(edge_extract_output$phylum2)
               d_shuffle<-with(edge_extract_output,table(phylum1,phylum2shuffle)) %>% as.data.frame()
               permute_vec[k]<-d_shuffle$Freq[which(d_shuffle$phylum1==d1[1,1] & d_shuffle$phylum2shuffle==d1[1,2])]
          } #loop end 
               print(paste("The most common co-ccurence was",d1[1,1],"&",d1[1,2]))
               
               twosidep<-mean(permute_vec<= min(d1[1,3],-d1[1,3])) + mean(permute_vec>= max(d1[1,3],-d1[1,3]))
               print(paste(" 2 tailed p value = ",twosidep))
               
               
         
  } #function end

2.5.2 Taxonomy of Bacteria-Fungi Interactions

#### MAMMALS
  mammal_bf<-edge_extract(ig.mammal)
  heatplot(mammal_bf) + scale_fill_gradient()

  occur_rand(mammal_bf,10000)
## [1] "The most common co-ccurence was Firmicutes & Ascomycota"
## [1] " 2 tailed p value =  0.7595"
#### BIRDS
  bird_bf<-edge_extract(ig.bird)
  
   heatplot(bird_bf) + scale_fill_gradient()

     occur_rand(bird_bf,10000)
## [1] "The most common co-ccurence was Actinobacteria & Ascomycota"
## [1] " 2 tailed p value =  0.0014"
#### INVERTS
  invert_bf<-edge_extract(ig.invert)
  heatplot(invert_bf) + scale_fill_gradient()

    occur_rand(invert_bf,10000)
## [1] "The most common co-ccurence was Proteobacteria & Ascomycota"
## [1] " 2 tailed p value =  0.6201"
#### FISH
    fish_bf<-edge_extract(ig.fish)
  heatplot(fish_bf) + scale_fill_gradient()

    occur_rand(fish_bf,10000)
## [1] "The most common co-ccurence was Proteobacteria & Basidiomycota"
## [1] " 2 tailed p value =  1"
### AMPHIBIA
  amphib_bf<-edge_extract(ig.amphib)
  heatplot(amphib_bf) + scale_fill_gradient()

    occur_rand(amphib_bf,10000)
## [1] "The most common co-ccurence was Proteobacteria & Ascomycota"
## [1] " 2 tailed p value =  1"
## Adjust P 
    p.adjust(c(0.765,0.002,0.6205,1,1))
## [1] 1.00 0.01 1.00 1.00 1.00

2.5.3 Combined Heatmaps for Composite Plot

mammal_heat<-heatplot(mammal_bf) + scale_fill_gradient() + labs(title="Mammalia",fill="Proportion") +
  theme(plot.title = element_text(hjust = 0.5,face="italic",size=15),legend.title = element_text(size=12))
bird_heat<-heatplot(bird_bf) + scale_fill_gradient()+ labs(title="Aves",fill="Proportion") +
  theme(plot.title = element_text(hjust = 0.5,face="italic",size=15),legend.title = element_text(size=12))
invert_heat<-heatplot(invert_bf) + scale_fill_gradient()+ labs(title="Insecta",fill="Proportion") +
  theme(plot.title = element_text(hjust = 0.5,face="italic",size=15),legend.title = element_text(size=12))
fish_heat<-heatplot(fish_bf) + scale_fill_gradient()+ labs(title="Actinopterygii",fill="Proportion") +
  theme(plot.title = element_text(hjust = 0.5,face="italic",size=15),legend.title = element_text(size=12))
amphibian_heat<-heatplot(amphib_bf) + scale_fill_gradient()+ labs(title="Amphibia",fill="Proportion") +
  theme(plot.title = element_text(hjust = 0.5,face="italic",size=15),legend.title = element_text(size=12))

heat_grid<-plot_grid(mammal_heat,bird_heat,invert_heat,fish_heat,amphibian_heat,nrow=3,ncol=2)
heat_grid

ggsave2('Heatmap Grid.pdf',heat_grid,width=10,height=10,units="in")

2.6 Species Level Network Stats from Co-Occur

#Read in Co-oCcur Data
    posneg<-read.csv('posneg.csv',header=T)
    head(posneg)
##                  Common_name Positive Negative
## 1 African palm weevil larvae       14        9
## 2                  Bank vole       43       25
## 3              Barred hamlet       21       24
## 4                   Blackcap      188      231
## 5            Capuchin monkey       78       64
## 6               Carrion crow       31       35
    posneg$Common_name<-as.character(posneg$Common_name)
    #posneg$Class<-metadata$Class[match(posneg$Common_name,metadata$Common_name)]
        metadata$Common_name<-as.character(metadata$Common_name)

#Combine With Metadata
    head(metadata)
##     Sample2             Species         Common_name      Taxa     Sex    Situ
## 1  Alytes-1 Alytes obstetricans Common midwife toad Amphibian Unknown Captive
## 2 Alytes-13 Alytes obstetricans Common midwife toad Amphibian Unknown Captive
## 3 Alytes-14 Alytes obstetricans Common midwife toad Amphibian Unknown Captive
## 4 Alytes-15 Alytes obstetricans Common midwife toad Amphibian Unknown Captive
## 5 Alytes-16 Alytes obstetricans Common midwife toad Amphibian Unknown Captive
## 6  Alytes-2 Alytes obstetricans Common midwife toad Amphibian Unknown Captive
##         Location Country CollectionYear  TissueStorage SampleType
## 1 London Zoo, UK Captive           2015 Frozen at -20C  Skin swab
## 2 London Zoo, UK Captive           2015 Frozen at -20C  Skin swab
## 3 London Zoo, UK Captive           2015 Frozen at -20C  Skin swab
## 4 London Zoo, UK Captive           2015 Frozen at -20C  Skin swab
## 5 London Zoo, UK Captive           2015 Frozen at -20C  Skin swab
## 6 London Zoo, UK Captive           2015 Frozen at -20C  Skin swab
##       ExtractionKit Environment    Class Order Rough_Class      Diet
## 1 Qiagen DNEasy kit  Freshwater Amphibia Anura    Amphibia Carnivore
## 2 Qiagen DNEasy kit  Freshwater Amphibia Anura    Amphibia Carnivore
## 3 Qiagen DNEasy kit  Freshwater Amphibia Anura    Amphibia Carnivore
## 4 Qiagen DNEasy kit  Freshwater Amphibia Anura    Amphibia Carnivore
## 5 Qiagen DNEasy kit  Freshwater Amphibia Anura    Amphibia Carnivore
## 6 Qiagen DNEasy kit  Freshwater Amphibia Anura    Amphibia Carnivore
    posnegj<- metadata %>% select("Common_name","Taxa","Class","Order")  %>% unique %>%  inner_join(x=posneg,"Common_name")
    posnegj<- mutate(posnegj,n_interactions=Positive + Negative)

#### MODELS - Binomial 
    posnegj_filter <- posnegj %>% filter(Class %in% c("Actinopterygii","Aves","Mammalia","Insecta","Amphibia"))

      degree_prop1<- brm(Positive|trials(n_interactions) ~ Class,family=binomial(),data=posnegj_filter)
## Compiling Stan program...
## Start sampling
      summary(degree_prop1)
##  Family: binomial 
##   Links: mu = logit 
## Formula: Positive | trials(n_interactions) ~ Class 
##    Data: posnegj_filter (Number of observations: 39) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept         0.46      0.17     0.11     0.80 1.01      774     1228
## ClassAmphibia    -0.20      0.20    -0.59     0.21 1.01      950     1789
## ClassAves        -0.13      0.18    -0.48     0.22 1.01      784     1189
## ClassInsecta     -0.60      0.23    -1.06    -0.16 1.00     1029     1843
## ClassMammalia     0.05      0.18    -0.30     0.39 1.01      813     1551
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
            pp_check(degree_prop1)
## Using 10 posterior samples for ppc type 'dens_overlay' by default.

                    conditional_effects(degree_prop1)
## Setting all 'trials' variables to 1 by default if not specified otherwise.

                          bayes_R2(degree_prop1)
##     Estimate   Est.Error      Q2.5     Q97.5
## R2 0.9314302 0.001351038 0.9284625 0.9336376

2.6.1 Model Plotting

## Taxa Names
  taxlabels<-c("Actinopterygii","Amphibia","Aves","Insecta","Mammalia")

########## Binomial Model 
    propmod_plotdata<-fixef(degree_prop1,summary=F)
    propmod_plotdata2<-data.frame(plogis(cbind(propmod_plotdata[,1],propmod_plotdata[,1]+propmod_plotdata[,2:5])))
    colnames(propmod_plotdata2)<-taxlabels
    
    propmod_plotdata3<- propmod_plotdata2 %>% pivot_longer(colnames(propmod_plotdata2),names_to="Class",values_to="value")
    propmod_plotdata3$facetval<-"proportion"
    
     ###Taxonomic Diffs
        # bird_mamm_prop<-propmod_plotdata2$Bird - propmod_plotdata2$Mammal
        #     mean(bird_mamm_prop); quantile(bird_mamm_prop,c(0.025,0.975))
        # mammal_amphib_prop<-propmod_plotdata2$Mammal - propmod_plotdata2$Amphibian
        #     mean(mammal_amphib_prop); quantile(mammal_amphib_prop,c(0.025,0.975))  
        #     
        #     
######### Plot - Proportion Positive Edges
     network_statplot2<- ggplot(propmod_plotdata3,aes(x=value,y=Class)) + geom_vline(xintercept = 0.5,lty="dashed") + theme_bw() + plotopts + stat_halfeye(shape=21,point_fill="white",point_size=5,colour="black",slab_alpha=0.5,aes(fill=Class)) + guides(fill=F)+ 
       labs(x="Proportion of Positive Edge Correlations")
     #ggsave2('Netowrk Edge Proportion Graph.pdf',network_statplot2,width=15,height=15,units="cm")
        #  ggsave2('Netowrk Edge Proportion Graph.jpeg',network_statplot2,width=15,height=15,units="cm")

2.6.2 Comparison of Proportion Positive vs Species Level Diversity

#Estimate Diversity - Fungi 
  network_rich_fungi<- myco %>% subset_samples(Common_name %in% posnegj$Common_name) %>% rarefy_even_depth(sample.size=500) %>% estimate_richness(measures = c("Shannon"))
## You set `rngseed` to FALSE. Make sure you've set & recorded
##  the random seed of your session for reproducibility.
## See `?set.seed`
## ...
## 51 samples removedbecause they contained fewer reads than `sample.size`.
## Up to first five removed samples are:
## Alytes-1Alytes-13Alytes-15Alytes-16Alytes-2
## ...
## 67OTUs were removed because they are no longer 
## present in any sample after random subsampling
## ...
network_rich_fungi$Common_name<- metadata$Common_name[match(rownames(network_rich_fungi),metadata$Sample2)]

#Estimate Diversity - Bactera 
  network_rich_bacteria<- micro %>% subset_samples(Common_name %in% posnegj$Common_name) %>% rarefy_even_depth(sample.size=500) %>% estimate_richness(measures = c("Shannon"))
## You set `rngseed` to FALSE. Make sure you've set & recorded
##  the random seed of your session for reproducibility.
## See `?set.seed`
## 
## ...
## 26 samples removedbecause they contained fewer reads than `sample.size`.
## Up to first five removed samples are:
## Alytes-15Alytes-16Bird25Bird29Bird62
## ...
## 17OTUs were removed because they are no longer 
## present in any sample after random subsampling
## ...
network_rich_bacteria$Common_name<- metadata$Common_name[match(rownames(network_rich_bacteria),metadata$Sample2)]

#Mean Values 
  network_rich_fungi_mean<- network_rich_fungi %>% group_by(Common_name) %>% summarise(mean_shannon=mean(Shannon))
  network_rich_bacteria_mean<- network_rich_bacteria %>% group_by(Common_name) %>% summarise(mean_shannon=mean(Shannon))

#Copy Across 
  posnegj$mean_fungal_shannon<-network_rich_fungi_mean$mean_shannon[match(posnegj$Common_name,network_rich_fungi_mean$Common_name)]
    posnegj$mean_bacteria_shannon<-network_rich_bacteria_mean$mean_shannon[match(posnegj$Common_name,network_rich_bacteria_mean$Common_name)]
  head(posnegj)
##                  Common_name Positive Negative         Taxa          Class
## 1 African palm weevil larvae       14        9 Invertebrate        Insecta
## 2                  Bank vole       43       25       Mammal       Mammalia
## 3              Barred hamlet       21       24         Fish Actinopterygii
## 4                   Blackcap      188      231         Bird           Aves
## 5            Capuchin monkey       78       64       Mammal       Mammalia
## 6               Carrion crow       31       35         Bird           Aves
##           Order n_interactions mean_fungal_shannon mean_bacteria_shannon
## 1    Coleoptera             23           0.3933196             2.3824486
## 2      Rodentia             68           0.6383061             2.4488759
## 3   Perciformes             45           0.7792424             1.8451177
## 4 Passeriformes            419           0.9932420             1.1894112
## 5      Primates            142           0.9083341             1.4398882
## 6 Passeriformes             66           0.6840271             0.6922338
  # Calculate Mean Proportion 
    posnegj$proportion<- with(posnegj,Positive / n_interactions)
  #Filter to The Big 5  
    posnegj_filter <- posnegj %>% filter(Class %in% c("Actinopterygii","Aves","Mammalia","Insecta","Amphibia"))
    
######### PLOTS - SPECIES LEVEL 
  proportion_shannon_bacteria_plot1<- ggplot(posnegj_filter,aes(x=mean_bacteria_shannon,y=proportion)) + geom_point(aes(color=Class))
    with(posnegj_filter,cor.test(mean_bacteria_shannon,proportion))
## 
##  Pearson's product-moment correlation
## 
## data:  mean_bacteria_shannon and proportion
## t = -0.35308, df = 37, p-value = 0.726
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3667596  0.2623658
## sample estimates:
##         cor 
## -0.05794821
  proportion_shannon_fungi_plot1<- ggplot(posnegj_filter,aes(x=mean_fungal_shannon,y=proportion)) + geom_point()
        with(posnegj_filter,cor.test(mean_fungal_shannon,proportion))
## 
##  Pearson's product-moment correlation
## 
## data:  mean_fungal_shannon and proportion
## t = 0.31016, df = 37, p-value = 0.7582
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.2689142  0.3606460
## sample estimates:
##        cor 
## 0.05092353
######### PLOTS - CLASS LEVEL 
    posneg_class<- posnegj_filter %>% group_by(Class) %>% summarise(mean_prop=mean(proportion),mean_bac_shannon=mean(mean_bacteria_shannon),mean_fungal_shannon=mean(mean_fungal_shannon))
    
  proportion_shannon_bacteria_classplot1<- ggplot(posneg_class,aes(x=mean_bac_shannon,y=mean_prop)) + geom_smooth(method="lm") + geom_point(shape=21,size=4,aes(fill=Class)) + labs(x="Mean Bacterial Shannon Diversity",y="Mean Proportion Positive Edges") + plotopts
  proportion_shannon_bacteria_classplot1
## `geom_smooth()` using formula 'y ~ x'

    with(posneg_class,cor.test(mean_bac_shannon,mean_prop))
## 
##  Pearson's product-moment correlation
## 
## data:  mean_bac_shannon and mean_prop
## t = 0.2051, df = 3, p-value = 0.8506
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8531904  0.9058763
## sample estimates:
##       cor 
## 0.1175935
  proportion_shannon_fungi_classplot1<- ggplot(posneg_class,aes(x=mean_fungal_shannon,y=mean_prop)) + geom_smooth(method="lm") + geom_point(shape=21,size=4,aes(fill=Class))+ labs(x="Mean Fungal Shannon Diversity",y="Mean Proportion Positive Edges") + plotopts
  proportion_shannon_fungi_classplot1
## `geom_smooth()` using formula 'y ~ x'

    with(posneg_class,cor.test(mean_fungal_shannon,mean_prop))
## 
##  Pearson's product-moment correlation
## 
## data:  mean_fungal_shannon and mean_prop
## t = 2.1157, df = 3, p-value = 0.1247
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3418874  0.9841717
## sample estimates:
##       cor 
## 0.7737782

2.7 Betweenness Randomisations

Here we randomly reassign node labels across kingdoms (bacteria & fungi) to test whether fungi are critical betweenness components

######### RANDOMISATIONS

  #Number of Shuffles
    nrand<-1000
            
###########MAMMALS 

  ###Betweenness of Fungi   
    mammal_between<-betweenness(ig.mammal)
    mammal_fungal_permute<- replicate(10000,mean(mammal_between[grep("fungi",sample(names(mammal_between)))]))
      mammal_b_true<-  mean(mammal_between[grep("fungi",names(mammal_between))])
          mean(mammal_fungal_permute>=mammal_b_true)*2
## [1] 0.6292
  #Random Removal of Bacteria or Fungi   
        mammal_nfungi<-length(which(tax_table(mammals)[,1]=="k__Fungi"))
        sample(seq(nrow(tax_table(mammals))),mammal_nfungi)
##  [1] 100 128 215  90  38  74 217  11  29  34 117 133  25 226 184 172 141  77 224
## [20]  20   4 126 156 169 137 183  44 101  58 132  10 160
        mammal_no_fungi<- delete_vertices(ig.mammal,which(tax_table(mammals)[,1]=="k__Fungi")) %>% cluster_fast_greedy() %>% modularity() 
         qr_mammal_fungidelete_perm <- replicate(10000, delete_vertices(ig.mammal,sample(seq(nrow(tax_table(mammals))),mammal_nfungi)),simplify = F) %>% lapply(cluster_fast_greedy) %>%
            sapply(modularity) 
         mean(mammal_no_fungi>=qr_mammal_fungidelete_perm) #0.0441
## [1] 0.0464
  ##Betweenness of Fungi       
      mammal_bf_bacbetween<- delete_vertices(ig.mammal,which(tax_table(mammals)[,1]=="k__Fungi")) %>% betweenness() %>% mean()
      qr_mammal_fungidelete_between_perm <- replicate(1000, delete_vertices(ig.mammal,sample(seq(nrow(tax_table(mammals))),mammal_nfungi)),simplify = F) %>% sapply(betweenness) %>% apply(2,mean)
         hist(qr_mammal_fungidelete_between_perm)

          mean(mammal_bf_bacbetween<=qr_mammal_fungidelete_between_perm)
## [1] 0.144
###########BIRDS 

  #Betweenness of Fungi           
    bird_between<-betweenness(ig.bird)
    bird_fungal_permute<- replicate(10000,mean(bird_between[grep("fungi",sample(names(bird_between)))]))
      bird_b_true<-  mean(bird_between[grep("fungi",names(bird_between))])
          mean(bird_fungal_permute<=bird_b_true)*2 #0.04
## [1] 0.047
    #Hist 
       hist(bird_fungal_permute)

2.7.1 Plot of Randomisations

#Assemble Data
  between_permute_data<-data.frame(randval=c(mammal_fungal_permute,bird_fungal_permute),Class=rep(c("Mammalia","Aves"),each=nrand))

# #Plots
  mammal_permplot1<- between_permute_data %>% filter(Class=="Mammalia") %>% ggplot(.,aes(randval)) + geom_density(fill="lightgray")
  mammal_permplot2<- mammal_permplot1<- mammal_permplot1 + theme_bw() + plotopts +  labs(x="Mammal Randomised Network Betweenness",y="Density")
  mammal_permplot3<- mammal_permplot2 + geom_vline(xintercept = mammal_b_true,cex=1.5,col="red")


  bird_permplot1<- between_permute_data %>% filter(Class=="Aves") %>% ggplot(.,aes(randval)) + geom_density(fill="lightgray")
  bird_permplot2<- bird_permplot1 + theme_bw() + plotopts +  labs(x="Bird Randomised Network Betweenness",y="Density")
  bird_permplot3<- bird_permplot2 + geom_vline(xintercept = bird_b_true,cex=1.5,col="red")

#Combine
  dummydf<-data.frame(Class=c("Mammalia","Aves"),realval=c(mammal_b_true,bird_b_true)) #real value 
  textdf<-data.frame(label= c("p = 0.044*","p = 0.63"),Class=c("Aves","Mammalia"),x=750,y=0.003)
  #textdf2<-data.frame(label= c("2-tailed p = 0.044","2-tailed p = 0.63"),Class=c("Aves","Mammalia"),x=850,y=0.003)

  perm_plot_all<-plot_grid(mammal_permplot3,bird_permplot3,labels=c("Mammalia","Aves"),ncol=2)
  permplot1<-ggplot(between_permute_data,aes(randval)) + geom_density(aes(fill=Class),alpha=0.5) + facet_wrap(.~Class)
  permplot2<- permplot1 + theme_bw() + plotopts  + guides(fill=F) + labs(x="Randomised Betweenness Value",y="Density")
  permplot3<- permplot2 + geom_vline(data=dummydf,aes(xintercept=realval),cex=1.5,color="red") +  geom_text(data = textdf,
            mapping = aes(x = x,
                          y = y,
                          label = label),size=8)

2.8 Bootstrap Sampling - Modularity Distributions

### Function That Strips Out  Fungal Nodes for Betweenness Calculations  
           fungus_subset<-function(matrix_in,physeq_obj){
              return(matrix_in[tax_table(physeq_obj)[,1]=="k__Fungi",])
           }  

### Function That Calculates True  Fungal Betweenness
          extract_fungal_between<-function(igraph_obj,physeq_obj){
            return(betweenness(igraph_obj)[tax_table(physeq_obj)[,1]=="k__Fungi"])
          }
           
############# Modularity and Group Calcs
    modularity_data<-data.frame(class=c("Mammalia","Aves","Insecta","Actinopterygii","Amphibia"))
    
  #Cluster Each IGraph PObject  
      mammal_cluster<-ig.mammal %>% cluster_fast_greedy() 
      bird_cluster<-ig.bird %>% cluster_fast_greedy() 
      insect_cluster<-ig.invert %>% cluster_fast_greedy() 
      fish_cluster<-ig.fish %>% cluster_fast_greedy() 
      amphib_cluster<-ig.amphib %>% cluster_fast_greedy() 
      
  #Strip Out Data
      modularity_temp<-c(modularity(mammal_cluster),modularity(bird_cluster),modularity(insect_cluster),modularity(fish_cluster),modularity(amphib_cluster))
      groups_temp<-c(max(mammal_cluster$membership),max(bird_cluster$membership),max(insect_cluster$membership),max(fish_cluster$membership),max(amphib_cluster$membership))
      components_temp<-c(components(ig.mammal)$no,components(ig.bird)$no,components(ig.invert)$no,components(ig.fish)$no,components(ig.amphib)$no)
      
    #Add TO Data Frame  
      modularity_data$modularity<-round(modularity_temp,3)
      modularity_data$groups<-groups_temp
      modularity_data$components<-components_temp
      
    #Stats
      with(modularity_data,cor(modularity,groups))
## [1] 0.3947702
            with(modularity_data,cor(modularity,components))
## [1] 0.6737694
    #Plot
      mod_groups_plot<-ggplot(modularity_data,aes(x=groups,y=modularity)) + geom_smooth(method="lm") + geom_point(size=5,shape=21,fill="white") + theme_bw() + plotopts + labs(x="Number of Groups",y="Modularity")
      
      mod_components_plot<-ggplot(modularity_data,aes(x=components,y=modularity)) + geom_smooth(method="lm") + geom_point(size=5,shape=21,fill="white") + theme_bw() + plotopts + labs(x="Number of Components",y="Modularity")
      
      mod_combined_plot<-plot_grid(mod_groups_plot,mod_components_plot,labels=c("A","B"),label_size=30)
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
        ggsave2('Modularity Combined Plot.pdf',mod_combined_plot,width=35,height=15,units="cm")
        mod_combined_plot

###########Modularity Distributions
  mammal_modularity_dist<-replicate(10000, delete_vertices(ig.mammal,sample(seq(nrow(tax_table(mammals))),round(length(V(ig.mammal))*0.1))),simplify = F) %>% lapply(cluster_fast_greedy) %>%
            sapply(modularity) 
 bird_modularity_dist<-replicate(10000, delete_vertices(ig.bird,sample(seq(nrow(tax_table(birds))),round(length(V(ig.bird))*0.1))),simplify = F) %>% lapply(cluster_fast_greedy) %>%
            sapply(modularity) 
  invert_modularity_dist<-replicate(10000, delete_vertices(ig.invert,sample(seq(nrow(tax_table(inverts))),round(length(V(ig.invert))*0.1))),simplify = F) %>% lapply(cluster_fast_greedy) %>%
            sapply(modularity) 
  fish_modularity_dist<-replicate(10000, delete_vertices(ig.fish,sample(seq(nrow(tax_table(fish))),round(length(V(ig.fish))*0.1))),simplify = F) %>% lapply(cluster_fast_greedy) %>%
            sapply(modularity) 
  amphib_modularity_dist<-replicate(10000, delete_vertices(ig.amphib,sample(seq(nrow(tax_table(amphibians))),round(length(V(ig.amphib))*0.1))),simplify = F) %>% lapply(cluster_fast_greedy) %>%
            sapply(modularity) 
  
  
######## Assemble Modularity Data
  modularity_plotdata<-data.frame(modularity=c(mammal_modularity_dist,bird_modularity_dist,invert_modularity_dist,fish_modularity_dist,amphib_modularity_dist),class=rep(c("Mammalia","Aves","Insecta","Actinopterygii","Amphibia"),each=10000))
           
  
######## Plot
  ggplot(modularity_plotdata,aes(x=modularity,y=class)) + stat_halfeye(shape=21,point_fill="white",point_size=5,colour="black",slab_alpha=0.5,aes(fill=class))

2.8.1 Combined Plot of Modularity, Randomisations and Proportion Positive Edges

######### Plot - Modularity 
    network_statplot1<-  ggplot(modularity_plotdata,aes(x=modularity,y=class)) + theme_bw() + plotopts + stat_halfeye(shape=21,point_fill="white",point_size=5,colour="black",slab_alpha=0.5,aes(fill=class)) + guides(fill=F)+ labs(x="Modularity",y="Class")
    
######### Plot - Proportion Positive Edges
     network_statplot2<- ggplot(propmod_plotdata3,aes(x=value,y=Class)) + geom_vline(xintercept = 0.5,lty="dashed") + theme_bw() + plotopts + stat_halfeye(shape=21,point_fill="white",point_size=5,colour="black",slab_alpha=0.5,aes(fill=Class)) + guides(fill=F)+ 
       labs(x="Proportion of Positive Edge Correlations")
     #ggsave2('Netowrk Edge Proportion Graph.pdf',network_statplot2,width=15,height=15,units="cm")
        #  ggsave2('Netowrk Edge Proportion Graph.jpeg',network_statplot2,width=15,height=15,units="cm")

    
    #Extract Colour Data fro Classes from Proportion Neg Plot
     g <- ggplot_build(network_statplot2)
      class_fill_cols<-unique(g$data[[2]]["fill"])[,1]
     
    #Rejig Permplot with custom fill
      classcols<-c("Mammalia"=class_fill_cols[5],"Aves"=class_fill_cols[3])
      permplot4<- permplot3 + scale_fill_manual(values=classcols)

2.8.2 Combined Plot With Modularity Distributions

## Part 1 - Subplot of Modularity, Betweeness Randomisations and Proportion Positive Edges 
     network_plot_sub1<-plot_grid(network_statplot2,network_statplot1,permplot4,labels=c("A","B","C"),label_size = 30,ncol=1,nrow=3)

#Part 2 - Add Heatmap 
  network_plot_sub2<- plot_grid(network_plot_sub1,heat_grid,ncol=2,labels=c("","D"),label_size=30,rel_widths = c(1.4,2))
  network_plot_sub2

#Save
   ggsave2('Combined Network Statistics Graph v4 Aug21.pdf',network_plot_sub2,width=42,height=30,units="cm")
      ggsave2('Combined Network Statistics Graph v4 Jul21.jpeg',network_plot_sub2,width=18,height=12,units="in")
     ggsave2('Combined Network Statistics Graph v4 Jul21.tiff',network_plot_sub2,width=18,height=12,units="in")